logo

Explore glacier_mass_change#

Initialize Python environment#

# Import standard libraries
import zipfile

# Import third party libraries
import cartopy.crs
import cdsapi
import IPython
import matplotlib.animation
import matplotlib.colors
import matplotlib.patches
import matplotlib.pyplot
import pandas
import xarray

# Ignore distracting warnings
import warnings
warnings.filterwarnings('ignore')

Download data#

We will use the Climate Data Store (CDS) API to download the glacier mass change dataset.

NOTE: To use the CDS API, you first need to register (if not already), find your UID and API key on your acount page, and fill them in below.

# Fill in your UID and API key (separated by a colon :)
KEY = 'UID:API Key'

c = cdsapi.Client(key=KEY)
c.retrieve(
    name='derived-gridded-glacier-mass-change',
    request={
        'variable': 'glacier_mass_change',
        'product_version': 'wgms_fog_2022_09',
        'format': 'zip',
        'hydrological_year': [
            '1975_76', '1976_77', '1977_78',
            '1978_79', '1979_80', '1980_81',
            '1981_82', '1982_83', '1983_84',
            '1984_85', '1985_86', '1986_87',
            '1987_88', '1988_89', '1989_90',
            '1990_91', '1991_92', '1992_93',
            '1993_94', '1994_95', '1995_96',
            '1996_97', '1997_98', '1998_99',
            '1999_20', '2000_01', '2001_02',
            '2002_03', '2003_04', '2004_05',
            '2005_06', '2006_07', '2007_08',
            '2008_09', '2009_10', '2010_11',
            '2011_12', '2012_13', '2013_14',
            '2014_15', '2015_16', '2016_17',
            '2017_18', '2018_19', '2019_20',
            '2020_21',
        ],
    },
    target='glacier_mass_change.zip'
)
2023-10-31 17:11:28,108 INFO Welcome to the CDS
2023-10-31 17:11:28,109 INFO Sending request to https://cds.climate.copernicus.eu/api/v2/resources/derived-gridded-glacier-mass-change
2023-10-31 17:11:28,507 INFO Request is completed
2023-10-31 17:11:28,509 INFO Downloading https://download-0012-clone.copernicus-climate.eu/cache-compute-0012/cache/data8/dataset-derived-gridded-glacier-mass-change-bca4231f-7fa2-4103-8bda-0165b3bd4f2d.zip to glacier_mass_change.zip (4.2M)
2023-10-31 17:11:46,985 INFO Download rate 233.2K/s 
Result(content_length=4410916,content_type=application/zip,location=https://download-0012-clone.copernicus-climate.eu/cache-compute-0012/cache/data8/dataset-derived-gridded-glacier-mass-change-bca4231f-7fa2-4103-8bda-0165b3bd4f2d.zip)

Since the data is downloaded as a zip file, we have to first unzip it.

with zipfile.ZipFile('glacier_mass_change.zip', 'r') as file:
    file.extractall('glacier_mass_change')

Read and inspect data#

The data are formatted as multiple netCDF4 files, one for each year, but they can be read as a single dataset using xarray.

ds = xarray.open_mfdataset('glacier_mass_change/*.nc4')
ds
<xarray.Dataset>
Dimensions:      (time: 46, lat: 360, lon: 720)
Coordinates:
  * time         (time) datetime64[ns] 1976-01-01 1977-01-01 ... 2021-01-01
  * lat          (lat) float64 89.75 89.25 88.75 88.25 ... -88.75 -89.25 -89.75
  * lon          (lon) float64 -179.8 -179.2 -178.8 -178.2 ... 178.8 179.2 179.8
Data variables:
    Glacier      (time, lat, lon) float64 dask.array<chunksize=(1, 360, 720), meta=np.ndarray>
    Uncertainty  (time, lat, lon) float64 dask.array<chunksize=(1, 360, 720), meta=np.ndarray>
Attributes:
    title:         Global gridded annual glacier mass change product
    project:       Copernicus Climate Change Service (C3S) Essential Climate ...
    data_version:  version-wgms-fog-2022-09
    institution:   World Glacier Monitoring Service - Geography Department - ...
    created_by:    Dr. Ines Dussaillant   ines.dussaillant@geo.uzh.ch
    references:    Fluctuation of Glagiers (FoG) database version wgms-fog-20...
    citation:      
    Conventions:   CF Version CF-1.8
    comment:       Brief data description:Horizontal resolution:\t0.5° (latit...

Each time value contains a date, but represents a hydrological year. For example, 1976-01-01 represents the 1975–1976 hydrological year. To simplify what follows, we reduce each time to the end year of the hydrological year.

ds['time'] = [date.year for date in ds['time'].values.astype('datetime64[D]').tolist()]
ds['time']
<xarray.DataArray 'time' (time: 46)>
array([1976, 1977, 1978, 1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987,
       1988, 1989, 1990, 1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999,
       2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011,
       2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021])
Coordinates:
  * time     (time) int64 1976 1977 1978 1979 1980 ... 2017 2018 2019 2020 2021

Plot mass change on a map#

Static map for one year#

We use matplotlib and cartopy to make a map of glacier mass change for one year.

# Select year
YEAR = 2021

# Create a map with a plate carrée projection
figure = matplotlib.pyplot.figure(figsize=(12, 6))
axis = matplotlib.pyplot.axes(projection=cartopy.crs.PlateCarree())

# Add title
axis.set_title(f"Hydrological year {YEAR - 1}{YEAR}")

# Add latitude, longitude gridlines
axis.gridlines(draw_labels=True, alpha=0.25, linestyle=':', color='white')

# Add white continents against a gray background
axis.set_facecolor('darkgray')
axis.add_feature(cartopy.feature.LAND, facecolor='lightgray')

# Define range and ticks of color scale
ticks = [-1, -0.1, -0.01, 0, 0.01, 0.1, 1]

# Plot mass loss in red and gain in blue
# A log scale is used for colors to highlight differences between regions
im = matplotlib.pyplot.pcolormesh(
    ds['lon'],
    ds['lat'],
    ds['Glacier'].loc[YEAR],
    cmap='RdBu',
    norm=matplotlib.colors.SymLogNorm(
        vmin=ticks[0], vmax=ticks[-1], linthresh=0.01
    )
)

# Add a colorbar with custom ticks
cbar = matplotlib.pyplot.colorbar(
    im,
    fraction=0.025,
    pad=0.05,
    label='Glacier mass change (Gt)',
    ticks=ticks
)
cbar.ax.set_yticklabels([f'< {ticks[0]}', *ticks[1:-1], f'> {ticks[-1]}'])

# Show plot
figure.tight_layout(pad=0)
matplotlib.pyplot.show()
../../_images/bf7355bf98e70b867859f1cd9a75c1342fbde034cb6ded2e126e3c18c128388d.png

Animation for all years#

Using the map above as a template, we create a map for each year and display them in an animation. Make sure to run this cell immediately after the cell above.

def animate(time_index: int) -> None:
    im.set_array(ds['Glacier'][time_index].values.ravel())
    year = ds['time'].values[time_index]
    axis.set_title(f"Hydrological year {year - 1}{year}")

animation = matplotlib.animation.FuncAnimation(
    figure,
    func=animate,
    frames=ds['time'].size,
    interval=500  # ms
)

IPython.display.HTML(animation.to_jshtml())

Plot timeseries of global mass change#

Compute mean of global change#

The global mass change is the sum of the mass changes reported for each spatial grid cell. We compute it by summing mass change across the latitude (lat) and longitude (lon) dimensions.

global_change_mean = ds.sum(dim=['lat', 'lon'])['Glacier']

Compute standard deviation of global change#

If we assume that the mass changes for which we computed the sum above are uncorrelated, the standard deviation of their sum is the square root of the sum of the squares of their standard deviatons:

\( f = \textstyle\sum_i m_i \\ \sigma_f = \sqrt {\textstyle\sum_i \sigma_i^2} \)

See Wikipedia: Propagation of uncertainty for more information.

We divide the uncertainty by 1.96 below because the authors report that this value is 1.96 \(\sigma\) (equivalent to a 95% confidence interval for a normal distribution).

global_change_std = ((ds['Uncertainty'] / 1.96) ** 2).sum(dim=['lat', 'lon']) ** 0.5

View timeseries as a table#

We use pandas to tabulate the mean and standard deviation computed above.

pandas.DataFrame({
    'year': ds['time'],
    'mean': global_change_mean,
    'std': global_change_std
})
year mean std
0 1976 90.711681 29.427761
1 1977 -9.916545 29.827922
2 1978 -7.919292 29.135105
3 1979 -20.050397 29.187522
4 1980 80.726975 29.364131
5 1981 86.211509 31.511226
6 1982 -26.793842 28.662013
7 1983 140.430476 28.653894
8 1984 -29.203053 28.564248
9 1985 32.698277 28.637171
10 1986 69.848798 27.594141
11 1987 134.985990 27.841201
12 1988 12.617256 27.050033
13 1989 -18.082409 26.847752
14 1990 -248.443006 27.656334
15 1991 -143.120973 26.829816
16 1992 133.672195 26.663257
17 1993 -87.847944 26.652052
18 1994 -89.210215 26.410071
19 1995 -200.446107 26.171569
20 1996 -83.870737 25.528025
21 1997 -327.845908 24.748657
22 1998 -158.643409 24.510523
23 1999 -229.454042 24.563616
24 2000 -144.428798 12.284889
25 2001 -162.298628 12.494367
26 2002 -160.765585 9.332476
27 2003 -190.294174 14.942672
28 2004 -379.090233 11.315378
29 2005 -403.113326 11.086587
30 2006 -285.326933 11.548195
31 2007 -340.458477 11.253159
32 2008 -174.674641 11.426205
33 2009 -257.799645 8.735710
34 2010 -297.322767 11.295674
35 2011 -390.143004 9.291715
36 2012 -231.262888 9.896819
37 2013 -254.336647 11.493179
38 2014 -165.881781 8.630415
39 2015 -219.722406 8.452308
40 2016 -378.749650 9.953543
41 2017 -273.700155 7.438282
42 2018 -217.651581 12.309113
43 2019 -500.817056 8.147666
44 2020 -398.682773 7.292266
45 2021 -336.283551 8.161970

Plot timeseries#

We use matplotlib to plot the mean and uncertainty of the global mass change for each year.

# Configure the figure
figure, axis = matplotlib.pyplot.subplots(1, 1, figsize=(9, 6))
axis.set_ylabel('Annual mass change (Gt)', fontsize=12)
axis.set_xlabel('Year', fontsize=12)

# Plot a horizontal line at 0 change
axis.axhline(y=0, alpha=0.5, linestyle=':', color='gray')

# Plot the uncertainty as 1.96 standard deviations below and above the mean
axis.fill_between(
    x=ds['time'],
    y1=global_change_mean - global_change_std * 1.96,
    y2=global_change_mean + global_change_std * 1.96,
    color='lightgray'
)

# Plot the mean as a line
axis.plot(
    ds['time'],
    global_change_mean,
    color='darkgray'
)

# Show figure
matplotlib.pyplot.show()
../../_images/351697d11a63264c6d92075da84819fe01cd62ed0c8fd69c07feb06983e6137b.png

We can also plot this in the style of the popular climate warming stripes.

# Configure the figure
figure, axis = matplotlib.pyplot.subplots(1, 1, figsize=(10, 1))
axis.set_xlim(0, global_change_mean.size)
axis.set_axis_off()

# Define a red (negative) to blue (positive) colorscale
colormap = matplotlib.pyplot.get_cmap('RdBu')

# Saturate colors at -500 and +500 Gt
normalizer = matplotlib.colors.Normalize(vmin=-500, vmax=500)

# Create a colored rectangle for each year
for i, value in enumerate(global_change_mean):
    axis.add_patch(matplotlib.patches.Rectangle(
        xy=(i, 0), width=1, height=1, facecolor=colormap(normalizer(value))
    ))

# Show figure
matplotlib.pyplot.show()
../../_images/54cd68b84f49836605a383083e8b92a52f46c9c871977a5c4e2506312659442f.png

Plot total global mass change#

Rather than plot the annual mass change, we would like to plot the total mass change accumulated from the beginning of the series to each subsequent year. To do so, we take the cumulative sum of the annual mass changes calculated above. To calculate the standard deviation of this total, we apply the same principle as before.

global_change_total = global_change_mean.cumsum()
global_change_total_std = (global_change_std ** 2).cumsum() ** 0.5

Again, we can use pandas to tabulate the total and its standard deviation.

pandas.DataFrame({
    'year': ds['time'],
    'total': global_change_total,
    'total_std': global_change_total_std
})
year total total_std
0 1976 90.711681 29.427761
1 1977 80.795136 41.901051
2 1978 72.875844 51.034816
3 1979 52.825448 58.791699
4 1980 133.552423 65.716939
5 1981 219.763931 72.881228
6 1982 192.970089 78.314650
7 1983 333.400565 83.392026
8 1984 304.197512 88.148433
9 1985 336.895789 92.683514
10 1986 406.744587 96.704035
11 1987 541.730577 100.632017
12 1988 554.347833 104.204161
13 1989 536.265424 107.607198
14 1990 287.822419 111.104373
15 1991 144.701446 114.297947
16 1992 278.373640 117.366733
17 1993 190.525696 120.354817
18 1994 101.315481 123.218399
19 1995 -99.130625 125.967158
20 1996 -183.001362 128.527837
21 1997 -510.847270 130.888888
22 1998 -669.490679 133.164059
23 1999 -898.944721 135.410627
24 2000 -1043.373519 135.966747
25 2001 -1205.672147 136.539612
26 2002 -1366.437732 136.858177
27 2003 -1556.731905 137.671508
28 2004 -1935.822139 138.135737
29 2005 -2338.935464 138.579920
30 2006 -2624.262397 139.060257
31 2007 -2964.720875 139.514833
32 2008 -3139.395516 139.981952
33 2009 -3397.195161 140.254267
34 2010 -3694.517928 140.708393
35 2011 -4084.660932 141.014849
36 2012 -4315.923820 141.361716
37 2013 -4570.260467 141.828164
38 2014 -4736.142249 142.090506
39 2015 -4955.864655 142.341679
40 2016 -5334.614305 142.689266
41 2017 -5608.314460 142.883010
42 2018 -5825.966041 143.412234
43 2019 -6326.783097 143.643494
44 2020 -6725.465870 143.828476
45 2021 -7061.749421 144.059877

Finally, we use matplotlib to plot the total and its uncertainty.

# Configure the figure
figure, axis = matplotlib.pyplot.subplots(1, 1, figsize=(9, 6))
start_year = ds['time'].values[0] - 1
axis.set_ylabel(f"Total mass change (Gt) since {start_year}", fontsize=12)
axis.set_xlabel('Year', fontsize=12)

# Plot a horizontal line at 0 change
axis.axhline(y=0, alpha=0.5, linestyle=':', color='gray')

# Plot the uncertainty as 1.96 standard deviations
axis.fill_between(
    x=ds['time'],
    y1=global_change_total - global_change_total_std * 1.96,
    y2=global_change_total + global_change_total_std * 1.96,
    color='lightgray'
)

# Plot the mean as a line
axis.plot(
    ds['time'],
    global_change_total,
    color='darkgray'
)

# Show figure
matplotlib.pyplot.show()
../../_images/3115556b90a60771ea6e0dfa51e91dc04f1d17080ce2d3e49d39f3d92df25abb.png